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Preface 


Energy  detection  in  the  RF  spectrum  is  the  most  basic  technique  for  signal 
detection.  Typically,  this  requires  establishing  an  energy  detection  threshold  based 
on  a  noise-only  condition  (i.e.,  no  signal  present  in  the  RF  spectrum).  Initially,  the 
area  of  exploration  for  this  report  was  to  examine  the  potential  of  an  automatic 
energy  detection  thresholding  algorithm  based  on  the  RF  measurement.  It  would 
not  require  the  preliminary  RF  spectrum  noise-only  measurements  to  establish  the 
energy  detection  threshold. 

A  series  of  RF  spectrum  measurements  were  collected  at  the  US  Army  Research 
Faboratory  in  2013.  The  data  files  contained  different  segments  of  the  RF  spectrum 
with  various  resolution  bandwidths.  In  addition,  the  data  files  represented  RF 
spectral  conditions  with  and  without  RF  signals.  Some  of  these  data  files  were  used 
to  examine  a  technique  for  automatic  energy  detection  thresholding  determination. 

This  report  summarizes  the  development  of  the  kurtosis  histogram  excision 
algorithm.  RF  spectrum  of  noise  only  is  considered  to  be  normal  or  Gaussian 
distributed.  The  histogram  function  converts  the  RF  spectrum  into  a  probability 
density  function.  The  threshold  is  then  determined  through  an  iterative  elimination 
of  the  high-bin  histogram  values  until  the  kurtosis  value  has  been  reached.  This  is 
the  fifth  of  5  reports  that  detail  the  energy  detection  techniques  examined  with  the 
recorded  RF  spectrum  measurements.1^ 


1  Tom  K.  An  automated  energy  detection  algorithm  based  on  morphological  filter  processing  with  a  semi¬ 
disk  structure.  Adelphi  (MD):  Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8271. 

2  Tom  K.  An  automated  energy  detection  algorithm  based  on  morphological  filter  processing  with  a 
modified  watershed  transform.  Adelphi  (MD):  Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL- 
TR-8270. 

Tom  K.  An  automated  energy  detection  algorithm  based  on  morphological  and  statistical  processing 
techniques.  Adelphi  (MD):  Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8272. 

4  Tom  K.  An  automated  energy  detection  algorithm  based  on  consecutive  mean  excision.  Adelphi  (MD): 
Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8268. 
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1.  Introduction 


Energy  detection  is  the  simplest  method  of  detecting  a  signal  in  the  frequency 
spectrum.  When  doing  so,  a  comparison  is  made  between  the  frequency  spectral 
component  energy  and  a  detection  threshold  level.  Knowledge  of  the  frequency 
spectrum  is  usually  necessary  to  establish  this  detection  threshold  level.  Various 
methods  can  be  used  to  determine  the  detection  threshold  level,  such  as  establishing 
the  system  noise  statistics  offline  and  setting  the  threshold  for  a  given  probability 
of  detection  versus  probability  of  false  alarm. 

The  goal  of  this  study  was  to  develop  an  algorithm  that  can  analyze  the  frequency 
spectrum  and  establish  a  threshold  value  based  on  the  spectral  components.  The 
automated  processing  employs  techniques  from  statistical  analysis  and  fault 
detection  used  in  detecting  ball-bearing  faults.  When  monitoring  vibrational  data 
from  a  ball  bearing,  a  fault  indicator  is  typically  set  to  a  value  of  3.5.  When  the  RF 
spectrum  does  not  contain  a  signal,  it  is  considered  to  be  noise  only.  The  RF  noise 
spectrum  is  typically  assumed  to  have  a  Gaussian  distribution.  For  a  Gaussian 
distribution,  the  kurtosis  has  a  value  of  3 .  The  algorithm  developed  sorts  the  spectral 
data  set  into  a  histogram  representation.  The  highest  bin  value  is  deleted  from  the 
histogram  until  the  remaining  data  set  has  a  kurtosis  value  corresponding  to  the 
desired  value. 

2.  Data  Collection  and  Statistical  Summary 

The  local  RF  spectrum  was  measured  in  2013  on  the  rooftop  of  building  204  at  the 
US  Army  Research  Faboratory’s  (ARF’s)  Adelphi  location.  An  Agilent  N9342CN 
spectrum  analyzer  and  a  Discone  antenna  was  used  to  collect  RF  spectrum  data. 
This  spectrum  analyzer  was  operated  under  the  control  of  a  Fab  VIEW  software 
program  to  acquire  and  store  data  with  different  resolution  bandwidth  (RBW)  from 
1  kHz  to  1  MHz. 

These  data  files  represent  various  sizes  of  RF  spectral  coverage  from  10  MHz  up 
to  4  GHz.  The  number  of  data  files  for  each  spectral  band  varied.  The  larger  RBW 
files  were  acquired  over  seconds  of  data  acquisition  time  versus  the  small  RBW 
files.  Small  RBW  provides  fine  spectral  resolution,  but  it  impacts  data  acquisition 
time  for  spectral  coverage  and  data  size.  Depending  on  the  RBW  and  spectral 
coverage,  a  data  file  could  require  a  few  hours  of  acquisition. 

Table  1  summarizes  the  data  collection  measurement  files.  The  RF  spectral  bands 
covered  the  spectrum  from  10  MHz  to  4  GHz.  In  general,  various  spectral  bands 
were  measured  with  4  RBW  configuration.  Data  file  size  is  inversely  proportional 
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to  the  RBW.  Data  file  size  is  proportional  to  the  spectral  band  coverage.  The  data 
size  varied  from  approximately  1  KSample  to  4  MSample  data  points  per  file. 


Table  1  Summary  of  the  RF  spectrum  collection 


Spectral  band  coverage 

Number  of  RBW 

measurement 

RBW 

1. 1-1.6  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

2-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

3-4  GHz 

3 

1  MHz 

10  kHz,  100  kHz,  1  MHz 

4-6  GHz 

3 

10  kHz,  100  kHz,  1  MHz 

10  MHz-1  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-2  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz^l  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

100  MHz-1  GHz 

1 

1  MHz 

100  kHz 

3.  Statistical  Processing 

3.1  Statistical  Analysis 

Statistical  analysis  is  the  mathematical  science  dealing  with  the  analysis  or 
interpretation  of  data.  The  data  analyst  uses  a  few  straightforward  statistical 
techniques  as  a  means  of  summarizing  the  collected  data.  These  statistical 
techniques  are  under  the  area  of  descriptive  statistics,  which  is  a  methodology  to 
condense  the  data  in  quantitative  terms. 

In  commercial  prognostics  and  diagnostic  vibrational  monitoring  applications, 
statistical  techniques  that  are  mainly  used  for  alarm  purposes  in  industrial  plants 
are  the  statistical  moments  of  order  2,  3,  and  4.  The  probability  density  function 
(PDF)  of  the  vibrational  time  series  of  a  good  bearing  has  a  Gaussian  distribution 
(also  known  as  a  normal  distribution),  whereas  a  damaged  bearing  results  in  a  non- 
Gaussian  distribution  with  dominant  tails  because  of  a  relative  increase  in  the 
number  of  high  levels  of  acceleration.  These  techniques  can  be  applied  to  the  RF 
spectral  data  with  a  different  interpretation  of  the  results.1 
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3.1.1  Histogram 

A  histogram  provides  a  tool  to  characterize  the  amplitude  of  a  series  of  data  into  a 
visualization  profile  that  can  relate  to  the  statistical  features.  The  calculation  is  as 
follows: 


Let  d  be  the  number  of  divisions  needed  to  divide  the  range  into,  let  /ijwith  0  < 
i  <  d  be  the  columns  of  the  histogram,  then 


n 


I 

7  =  0 


~ri(xi)>  Vt,  0  <  i  <  d 
n 


nix)  = 


1 1  i/  i(jnax(xi)-min(xi ))  ^  (/+l)(max(xi)-rmn(xj)) 


d 


d 


0,  otherwise 


The  histogram  upper  bound  (/%)  and  lower  bound  ( hL )  are  defined  as 

A 

hu  =  max{xf)  + 


A 

hL  =  max(xi)  -  - , 


where  A  =  max(xj)  —  min 


(  xt 

V[n  -  1] 


■) 


3.1.2  Moments 


If  these  moments  are  calculated  about  the  mean,  they  are  called  central  statistical 
moments.  The  first  and  second  moments  are  well  known,  being  the  mean  and  the 
variance,  respectively.  These  are  analogous  to  the  first  and  second  area  moments  of 
inertia  with  the  area  shape  defined  by  the  PDF.  The  third  moment  is  termed 
skewness  and  the  fourth  moment  is  termed  kurtosis.  The  general  equation  for  the 
order  of  moment  is  as  follows: 


Mr 


N 

i= 1 


where  p  is  the  order  of  the  moment, 
N  is  the  number  of  data  value, 
i  is  the  index  of  the  data  value,  and 
x  is  the  mean  value  of  the  data  set. 
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3.1.3  Mean 

Mean  is  the  most  common  measure  of  a  statistical  distribution.  In  this  case,  mean 
is  the  arithmetic  average  for  a  set  of  measurements. 


i= 1 


3.1.4  Variance 

Variance  is  a  measure  of  the  dispersion  of  a  waveform  about  its  mean — also  called 
the  second  moment  of  the  measurements. 

N 

ff2 = " 502  ■ 

1= 1 


3.1.5  Standard  Deviation 

Standard  deviation  is  a  measure  of  the  variation  of  a  set  of  data  values.  The  standard 
deviation  is  defined  as  the  square  root  of  the  variance  moment. 


a 


\ 

1 

i= 1 


(x i  -  x)2 


3.1.6  Kurtosis 

Kurtosis  is  the  fourth  statistical  moment,  normalized  by  the  standard  deviation  to 
the  fourth  power.  It  is  a  measure  of  whether  the  data  are  peaked  or  flat  relative  to  a 
normal  distribution.  The  noise  in  the  RF  spectrum  is  typically  considered  to  have  a 
normal  distribution.  The  normal  distribution  has  a  value  of  3. 

M4 

K  =  — -  . 

er4 

-  *)4 

K~  a* 

N 

K  =  a r^'L(Xi  ~ s)4  ■ 

i= l 
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3.1.7  Maximum 


“Max”  is  the  largest  value  of  a  set  of  numbers. 

y  —  max[x(n)]. 


3.1.8  Minimum 

“Min”  is  the  smallest  value  of  a  set  of  numbers. 

y  —  min[x(n)]. 


3.1.9  Median 

The  statistical  median  is  an  order  statistic  that  gives  the  “middle”  value  of  a  set  of 
samples.  Median  is  the  middle  value  of  a  set  of  data  values  that  divides  the  set  into 
2  groups.  Half  the  groups  exist  below  and  half  exist  above  this  value. 

3.1.10  Crest  Factor  (CF) 

Crest  factor  (CF)  is  a  measure  of  a  waveform  showing  the  ratio  of  peak  values  to 
the  effective  value.  In  other  words,  CF  indicates  how  extreme  the  peaks  are  in  a 
waveform. 

Crest  Factor  —  ^  ^peak  _ 

^rms 

Noise  sources  are  characterized  by  their  CF,  which  is  the  peak  to  average  ratio  of 
the  noise.  In  a  technical  bulletin,  XiTRON  reported  CF  values  between  5  and  7  for 
random  noise.2  For  example,  a  5: 1  CF  of  the  noise  voltage  is  201og5  =  14  dB.  This 
is  a  measure  of  the  quality  of  the  noise  distributions  and  one  way  to  measure  its 
Gaussian  nature.3  For  the  purpose  of  algorithm  development,  the  CF  equation  was 
modified  as  follows: 

Max 

Crest  Factor  —  - - - . 

Median 

3.2  Statistical  Summary 

There  were  31  different  groupings  of  the  RF  spectrum  data  measurements.  The 
number  of  data  files  under  each  of  the  main  groupings  was  varied.  For  the  purpose 
of  developing  the  algorithm,  only  a  single  data  file  was  selected  from  each  group. 
Each  data  file  was  processed  to  obtain  the  following  characteristics:  RBW,  mean, 
standard  deviation,  median,  max,  min,  kurtosis,  and  CF.  A  summary  of  the  results 
is  in  Table  2. 
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The  following  results  were  noted: 

•  The  smaller  the  RBW,  the  lower  the  noise  floor.  The  equation  for  thermal 
noise  power  is  P  —  kTB,  where  B  is  bandwidth.  In  this  case,  the  RBW  of  1 
kHz  has  the  lowest  noise  value.  The  RBW  of  1  MHz  has  the  highest  noise 
value. 

•  Each  10-fold  increase  in  bandwidth  results  in  a  10-dB  increase  in  noise 
power.  This  relationship  is  illustrated  in  this  data  set. 

•  The  calculated  mean  and  median  values  are  very  close  for  a  given  RF 
measurement  configuration. 

•  A  visual  examination  of  the  RF  data  files  shows  that  there  are  2  categories 
of  data:  signal  plus  noise  and  noise  only.  The  kurtosis  value  correlates  with 
the  2  categories.  Typical  kurtosis  values  for  the  noise-only  files  were  3  to 
3.6. 


Table  2  Statistical  analysis  of  the  RF  spectrum  measurements 


Filename 

Spectral  Band 

RBW 

Mean 

SD 

Median 

Max 

Min 

Kurtosis 

CF 

A  i  r_te  st_l.  1G  H  z_l.  6G  H  z  b_03_28_14_06_33_22 

1.1-  1.6GHz 

1000 

-112.694 

3.8444 

-112.3 

-101.2 

-138.4 

3.5571 

11.1 

Ai  r_test_l.  lGHz_l.  6GHza_03_27_14_07_14_05 

1.1-  1.6GHz 

10000 

-102.12 

3.7226 

-101.8 

-91.69 

-124.5 

3.5134 

10.11 

Ai  r_test_l.  lGHz_l.  6GHzc_03_31_14_06_47_26 

1.1-  1.6GHz 

100000 

-90.9109 

2.597 

-90.75 

-83.13 

-101.2 

3.0217 

7.62 

Ai  r_test_l.lGHz_l.  6GHzd_04_01_14_06_54_03 

1.1-  1.6GHz 

1000000 

-80.5634 

3.5741 

-80.745 

-55.36 

-89.8 

14.6861 

25.385 

Ai  r_te  st_2G  H  z_3G  H  za_06_05_14_07_15_58 

2- 3 GHz 

1000 

-111.219 

4.4991 

-111.1 

-77.55 

-139 

6.8918 

33.55 

A  i  r_te  st_2G  H  z_3G  H  z  b_05_29_14_06_28_28 

2- 3 GHz 

10000 

-101.006 

4.1808 

-100.8 

-74.67 

-122 

6.1995 

26.13 

Ai  r_te  st_2G  H  z_3G  H  zc_05_29_14_04_09_27 

2- 3 GHz 

100000 

-89.5666 

3.4154 

-89.68 

-64.27 

-100.3 

13.2282 

25.41 

Ai  r_test_2G  Hz_3G  Hzd_05_28_14_00_06_18 

2- 3 GHz 

1000000 

-79.6014 

3.0753 

-79.6 

-58.74 

-88.88 

8.4466 

20.86 

A  i  r_te  st_3G  H  z_4G  H  z  b_06_12_14_06_30_17 

3-4  GHz 

10000 

-100.51 

3.7072 

-100.2 

-89.68 

-123.7 

3.5901 

10.52 

A  i  r_te  st_3G  H  z_4G  H  zC_06_l  1_14_06_44_40 

3-4  GHz 

100000 

-89.1625 

2.5513 

-89.075 

-80.68 

-102.3 

3.0774 

8.395 

A  i  r_te  st_3G  H  z_4G  H  zd_06_09_14_07_09_5 1 

3-4  GHz 

1000000 

-79.0828 

2.445 

-78.94 

-72.76 

-88.5 

2.9651 

6.18 

A  i  r_te  st_4G  H  z_6G  H  za_04_24_14_07_05_26 

4- 6 GHz 

1000 

-107.559 

4.1793 

-107.2 

-94.31 

-136.3 

3.3696 

12.89 

A  i  r_te  st_4G  H  z_6G  H  z  b_04_28_14_06_23_12 

4- 6 GHz 

10000 

-97.245 

4.0398 

-96.95 

-84.05 

-119.7 

3.3249 

12.9 

A  i  r_te  st_4G  H  z_6G  H  zc_04_29_14_06_48_42 

4- 6  GHz 

100000 

-85.9125 

3.082 

-85.77 

-75.84 

-98.78 

2.9068 

9.93 

A  i  r_te  st_10M  H  z_lG  Hza_04_17_14_07_09_30 

10  MHz  - 1  GHz 

1000 

-110.348 

8.6126 

-112.1 

-28.37 

-140.3 

6.4788 

83.73 

A  i  r_te  s  t_10M  H  z_lG  H  z  b_04_2 1_14_07_17_44 

10  MHz  - 1  GHz 

10000 

-100.168 

8.4775 

-102 

-27.95 

-126.7 

6.7602 

74.05 

A  i  r_te  st_10M  H  z_lG  H  zc_04_22_14_06_39_58 

10  MHz  - 1  GHz 

100000 

-87.948 

8.8113 

-90.74 

-27.28 

-103.6 

6.513 

63.46 

Air_test_10MHz_lGHzd_04_23_14_06_36_16 

10  MHz  - 1  GHz 

1000000 

-78.1054 

8.3327 

-80.59 

-27.82 

-89.88 

7.5282 

52.77 

A  i  r_te  st_10M  H  z_2G  H  zd_03_06_14_08_52_04 

10  MHz -2  GHz 

1000 

-111.275 

6.9798 

-112.1 

-27.34 

-138.3 

9.5085 

84.76 

A  i  r_te  st_10M  H  z_2G  H  zC_02_27_14_06_43_55 

10  MHz -2  GHz 

10000 

-100.894 

7.0157 

-101.8 

-26.58 

-125.9 

10.4224 

75.22 

A  i  r_te  st_10M  H  z_2G  H  z  B_02_25_14_06_35_53 

10  MHz -2  GHz 

100000 

-89.242 

6.7487 

-90.59 

-26.09 

-102.4 

12.7961 

64.5 

A  i  r_te  st_10M  Hz_2G  Hz_02_20_14_06_55_36 

10  MHz -2  GHz 

1000000 

-78.8448 

7.3203 

-80.6 

-26.54 

-92.41 

11.3819 

54.06 

Air_test_10MHz_3GHz_03_ll_14_06_32_51 

10  MHz -3  GHz 

1000 

-111.373 

6.1356 

-111.8 

-29.42 

-140.5 

11.0474 

82.38 

Ai  r_te  st_10MHz_3G  Hz  B_03_13_14_06_32_45 

10  MHz -3  GHz 

10000 

-100.787 

6.3591 

-101.4 

-25.36 

-125.9 

11.7254 

76.04 

A  i  r_t  e  s  t_10  M  H  z_3G  H  z  C_03_18_14_06_32_08 

10  MHz -3  GHz 

100000 

-89.1614 

6.0257 

-90.14 

-24.16 

-102.4 

15.3376 

65.98 

A  i  r_t  e  s  t_10  M  H  z_3G  H  z  d_03_20_14_08_09_07 

10  Mhz  -  3  GHz 

1000000 

-78.6492 

6.7418 

-80.8 

-25.81 

-90.38 

13.0673 

54.99 

Air_test_10MHz_4GHza_04_10_14_07_01_22 

10  MHz -4  GHz 

1000 

-111.158 

5.6813 

-111.4 

-28.22 

-143.1 

10.9573 

83.18 

Air_test_10MHz_4GHzb_04_ll_14_06_10_18 

10  MHz -4  GHz 

10000 

-100.838 

5.5851 

-101.1 

-26.85 

-128 

12.2576 

74.25 

Air_test_10MHz_4GHzc_04_15_14_07_01_09 

10  MHz -4  GHz 

100000 

-89.4431 

5.1076 

-90.03 

-25.05 

-103.5 

18.317 

64.98 

A  i  r_te  st_10M  Hz_4G  Hzd_04_16_14_06_42_27 

10  MHz -4  GHz 

1000000 

-79.0396 

5.5441 

-79.89 

-26.72 

-91.64 

16.9073 

53.17 

Ai  r_te  st_100M  H  z_lG  Hz_02_19_14_07_20_46 

100  MHz  - 1  GHz 

100000 

-88.0772 

9.0701 

-90.82 

-23.91 

-102.2 

6.994 

66.91 
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4.  Algorithm 


This  development  process  was  executed  in  MATLAB.  Appendix  A  is  the  code  used 
to  process  and  generate  the  enclosed  RF  spectrum  signature  with  the  corresponding 
results  of  the  detection  threshold  level.  Basically,  the  algorithm  generates  a 
histogram  of  the  RF  spectrum  file  and  examines  the  kurtosis  statistics  of  the 
histogram.  If  the  kurtosis  is  less  than  3.6,  the  threshold  is  declared  as  the  maximum 
value  of  the  signature  file.  If  it  is  greater  than  3.6,  the  histogram  bin  information  is 
sequentially  eliminated.  Starting  with  the  highest  bin  value,  the  bins  are  deleted 
from  the  histogram  until  the  kurtosis  of  the  reduced  histogram  is  less  than  3.6.  The 
following  is  a  description  of  the  detection  thresholding  level  generation  based  on 
this  technique: 

1)  Determine  some  statistical  parameters  of  the  spectral  data  (i.e.,  min,  max, 
and  kurtosis).  Use  these  parameters  (min  and  max)  to  determine  the 
dynamic  range  of  the  spectral  amplitude  data. 

2)  Generate  a  histogram  of  the  spectral  amplitude  data.  The  number  of  bins  for 
the  histogram  is  twice  the  dynamic  range,  as  determined  in  step  1. 

3)  Store  the  results  of  the  histogram  sort  into  2  arrays:  Amplitudecount  and 
Amplitudevalue. 

4)  Generate  a  new  array  called  Amplitudediff.  The  new  array  is  generated  by 
taking  the  derivative  of  the  Amplitudecount  array. 

5)  Determine  the  smallest  value  in  the  Amplitudediff  array.  Then  determine 
the  location  index  in  the  Amplitudediff  array  that  corresponds  to  this 
smallest  value. 

6)  Search  for  the  bin  that  corresponds  to  the  smallest  value  in  the 
Amplitudediff  array.  Establish  the  index  to  start  the  search  for  kurtosis  that 
matches  the  criteria.  Calculate  the  index  count  as  the  smallest 
Amplitudediff  index  +  [0.6  x  (max  index  -  smallest  Amplitudediff 
index)].  This  step  is  performed  to  reduce  the  threshold  search  time. 

7)  If  the  initial  kurtosis  value  is  less  than  3.6,  then  use  the  max  value  as  the 
threshold. 

8)  If  the  initial  kurtosis  value  is  greater  than  3.6,  form  a  shortened  copy  of  the 
histogram  array.  Using  the  Amplitude  value  array,  form  a  new  array  by 
deleting  the  values  starting  at  the  index  count  from  step  6. 
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9)  Calculate  the  kurtosis  on  the  new  array.  If  this  new  kurtosis  value  is 
greater  than  3.6,  repeat  step  8  by  reducing  the  index  count  by  1.  Repeat 
this  until  the  kurtosis  is  less  than  3.6.  The  bin  value  at  this  index  is  then 
selected  to  be  the  threshold  value. 

5.  Conclusion 


The  kurtosis-histogram  excision  algorithm  works  well.  Appendix  A  is  the 
MATLAB  code  developed  for  this  algorithm,  which  was  based  on  prior  work  on 
vibrational  analysis  of  ball  bearings.  In  vibrational  analysis,  a  typical  kurtosis  value 
of  approximately  3.5  is  used  as  a  threshold  to  indicate  that  there  is  a  deviation  from 
a  normal  distribution.  In  this  algorithm,  the  histogram  processing  reduces  the  RF 
spectrum  data  to  a  form  where  it  can  be  compared  to  a  normal  distribution.  In  this 
case,  a  normal  distribution  represents  a  noise-only  condition.  This  iterative 
algorithm  starts  by  eliminating  the  high  bin  values  in  the  histogram  array  and 
recalculating  the  kurtosis  value.  When  the  kurtosis  value  is  less  than  3.6,  then  the 
histogram  bin  value  is  declared  as  the  RF  spectrum  threshold. 

Appendix  B  displays  the  results  of  threshold  overlaid  on  top  of  the  RF  spectrum 
signature.  The  graphs  are  intended  to  provide  the  reader  with  a  qualitative  sense  of 
the  effectiveness  of  the  automatic  detection  threshold  generation  algorithm.  The 
algorithm  works  well  for  the  RF  spectrum  where  there  are  signals  present  as  well 
as  for  the  RF  spectrum  that  was  judged  to  be  noise  only.  A  visual  inspection  of  the 
graphs  shows  that  the  algorithm  works  well  in  minimizing  the  number  of  false 
alarms  while  permitting  the  detection  of  signals.  Examination  of  the  graphs  shows 
that  the  spectrum  is  not  necessarily  flat  over  the  spectral  bands.  This  algorithm 
results  in  a  single  value  threshold  detection  level  across  the  spectral  bands. 
Therefore,  in  some  cases,  the  detection  of  certain  signals  may  be  degraded  since  the 
detection  threshold  is  higher  than  necessary. 
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Appendix  A.  MATLAB  Code 
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function  KurtosisThreshold2report () 


%  Make  selection  for  data  file  to  process 
%  SELECT  =  1  to  31 

SELECT  =  31; 


dir{l}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHza\  '  ; 

filename{l}  =  ' Air_test_l . lGHz_l . 6GHza_03_27_14_07_14_05 ' ; 

dir {2}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzb\ '  ; 

filename{2}  =  ' Air_test_l . lGHz_l . 6GHzb_03_28_14_06_33_22 ' ; 

dir{3}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzc\ '  ; 

filename{3}  =  'Air_test_l . lGHz_l . 6GHzc_03_31_14_0 6_47_2 6 ' ; 

dir{4}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzd\  '  ; 

filename{4}  =  'Air_test_l . lGHz_l . 6GHzd_04_01_14_06_54_03 ' ; 

dir{5}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHza\ ' ; 

filename{5}  =  'Air_test_2GHz_3GHza_06_05_14_07_15_58 ' ; 

dir{6}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzb\ ' ; 

filename{6}  =  ' Air_test_2GHz_3GHzb_0 5_2 9_1 4_0 6_2 8_2 8 ' ; 

dir{7}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzc\ '  ; 

filename{7}  =  'Air_test_2GHz_3GHzc_05_29_14_04_09_27 ' ; 

dir{8}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzd\ ' ; 

filename{8}  =  'Air_test_2GHz_3GHzd_05_28_14_00_06_18 ' ; 

dir{9}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_3GHz_4GHzb\ ' ; 

filename{9}  =  'Air_test_3GHz_4GHzb_06_12_14_06_30_17 ' ; 
dir{10}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzC\ '  ; 
filename{10}  =  'Air_test_3GHz_4GHzC_06_ll_14_06_44_40 ' ; 

dir{ll}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzd\ ' ; 
filename{ll}  =  ' Air_test_3GHz_4GHzd_06_09_14_07_09_51 ' ; 


Approved  for  public  release;  distribution  is  unlimited. 


12 


dir{12}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHza\ ' ; 
filename{12}  =  ' Air_test_4GHz_6GHza_04_24_14_07_05_26 ' ; 

dir{13}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzb\ ' ; 
filename{13}  =  ' Air_test_4GHz_6GHzb_04_28_14_06_23_12 '  ; 

dir{14}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzc\ ' ; 
filename{14}  =  ' Air_test_4GHz_6GHzc_04_29_14_06_48_42 ' ; 

dir{15}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c085\Air_test_10MHz_lGHza\ ' 
filename{15}  =  ' Air_test_10MHz_lGHza_04_17_14_07_09_30 ' ; 

dir{16}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_lGHzb\ ' 
filename{16}  =  ' Air_test_10MHz_lGHzb_04_21_14_07_17_44 ' ; 

dir{17}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_lGHzc\ ' 
filename{17}  =  'Air_test_10MHz_lGHzc_04_22_14_06_39_58 ' ; 

dir{18}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_lGHzd\ ' 
filename{18}  =  ' Air_test_10MHz_lGHzd_04_23_14_06_36_16 ' ; 

dir{19}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHz\ ' ; 
filename{19}  =  'Air_test_10MHz_2GHz_02_20_14_06_55_36 '  ; 

dir{20}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHzB\ ' 
filename{20}  =  ' Air_test_10MHz_2GHzB_02_25_14_06_35_53 ' ; 

dir{21}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_2GHzC\ ' 
filename{21}  =  ' Air_test_10MHz_2GHzC_02_27_14_06_43_55 ' ; 

dir{22}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c085\Air_test_10MHz_2GHzd\ ' 
filename{22}  =  ' Air_test_10MHz_2GHzd_03_06_14_08_52_04 ' ; 

dir{23}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHz\ ' ; 
filename{23}  =  ' Air_test_10MHz_3GHz_03_ll_14_06_32_51 ' ; 

dir{24}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_3GHzB\ ' 
filename{24}  =  'Air_test_10MHz_3GHzB_03_13_14_06_32_45 ' ; 

dir{25}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzC\ ' 
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=  ' Air_test_10MHz_3GHzC_03_18_14_06_32_08 ' ; 
dir{26}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzd\ ' ; 
filename{26}  =  ' Air_test_10MHz_3GHzd_03_20_14_08_09_07 ' ; 

dir{27}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHza\ ' ; 
filename{27}  =  ' Air_test_10MHz_4GHza_04_10_14_07_01_22 ' ; 

dir{28}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzb\ ' ; 
filename{28}  =  ' Air_test_10MHz_4GHzb_04_ll_14_06_10_18 ' ; 

dir{29}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzc\ ' ; 
filename{29}  =  ' Air_test_10MHz_4GHzc_04_15_14_07_01_09 ' ; 

dir{30}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzd\ ' ; 
filename{30}  =  ' Air_test_10MHz_4GHzd_04_16_14_06_42_27 ' ; 

dir{31}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_100MHz_lGHz\ ' ; 
filename{31}  =  ' Air_test_100MHz_lGHz_02_19_14_07_20_46 ' ; 

str_meta  = 

sprint f ( ' %s%s . mspecrawdata ' , dir { SELECT } , filename { SELECT } ) 
str_data  = 

sprint f ( ' %s%s . specrawdata ' , dir { SELECT } , filename { SELECT } ) ; 

%  Get  metadata  on  selected  data  file 

fid__meta  =  f open  ( str_meta)  ; 

META  =  textscan ( f id_meta,  ' %s  '  ) ; 

ave  =  META{ 1 } {!}; 

ref  =  META{ 1 } {2 }  ; 

attn  =  META{ 1 } { 3 }  ; 

rbw  =  META{ 1 } {4}; 


%  Read  in  data  file 


fid_data  =  f open ( str_data) ; 

g=0; 

f=0; 

a=0 ; 

while (g==0 ) 

ft=fgetl (fid_data)  ; 
f = [ f , str2num ( ft )  ]  ; 


at=fgetl (fid_data)  ; 
a= [ a, str2num (at ) ] ; 
g=feof (fid_data)  ; 

end 


Approved  for  public  release;  distribution  is  unlimited. 


14 


f=f (2 : end) ; 
a=a (2 : end) ; 


%  Plot  data  file 

figure  ( 1 ) ; 
plot (f /le6, a) ; 
axis  tight; 
grid; 

xlabel ( 'Frequency  (MHz)  ' )  ; 

ylabel ( 'Power  (dBm) ' ) ; 

title (strcat (rbw,  '  ' , ave,  '  -  -  ' , ref,  '  - 

' , attn) , ' Interpreter ' , ’ none ' ) ; 

%  Determine  some  statistics  values  for  data 

M  =  mean (a) ; 

Med  =  median (a); 

S  =  std (a) ; 

Max  =  max  (a) ; 

Min  =  min  (a) ; 

Kurt  =  kurtosis  (a) ; 

Range  =  abs (Min  -  Max) ; 

Number  =  floor (Range ) ; 

Bins  =  Number  *  2; 

disp ( [M  S  Med  Max  Min  Kurt]); 


[Amplitudearray  X]  =  hist (a, Bins) ; 

%  Plot  sorted  amplitdude  data  array 
figure  (2 ) ; 

plot (Amplitudearray) ; 


%  This  algorithm  uses  the  histogram  and  kurtosis  processing  to 
determine 
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%  the  threshold.  The  RF  spectrum  is  sorted  by  the  histogram. 
There  is  an 

%  iterative  elimination  of  the  highest  values  from  the  histogram 
until  the 

%  criteria  of  the  kurtosis  of  the  newly  formed  array  is  less  than 
3.6.  The 

%  value  of  3.6  is  somewhat  an  empirical  obervation  from 
vibrational 

%  analysis  to  determine  when  the  signal  has  deviated  from  a 
normal 

%  distribution  or  noise  only  condition, 
index  =  Bins; 


[Amplitudecount  Amplitudevalue]  =  hist (a, Bins) ; 
Histovalue  =  a; 
indexvalue  =  1; 

%  Take  derivative  of  the  histogram  values  and  plot  array 
figure ( 4 ) ; 

Amplitudedif f  =  diff (Amplitudecount) ; 
plot  (Amplitudedif f)  ; 

%  Find  the  smallest  value  from  the  derivatives 
Amplitudedif fmin  =  min (Amplitudedif f)  ; 

%  Determine  starting  index  of  histogram  array  to  start 
elimination 

for  indexcount  =  1: index-1 

if  Amplitudedif f ( (indexcount) )  ==  Amplitudedif fmin 
indexsubtraction  =  index  -  floor ((index  - 
indexcount) * . 6) ; 
end 


end 


for  m  =  1 : index 

for  n=  1 : Amplitudecount (m)  ; 

Histovalue (indexvalue)  =  Amplitudevalue (m) ; 
indexvalue  =  indexvalue  +  1; 

end 

end 

Histolength  =  length (Histovalue)  ; 

ProdKurt  =  kurtosis (Histovalue) ; 
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%  Set  the  starting  point  for  elimination  from  the  histogram 
array 

index  =  indexsubtraction; 

%  Check  to  see  if  kurtosis  value  of  the  newly  formed  histogram  is 
below 

%  3.6.  It  criteria  has  been  meet,  then  the  threshold  is  the 
value  of  the 
%  histogram  bin. 

if  ProdKurt  <  3.6 

index  =  Bins; 
else 

%  While  kurtosis  value  is  greater  than  3.6,  eliminate  next  bin 
value  and 

%  form  reduced  array  for  kurtosis  calculation 

while  ProdKurt  >  3.6 

Histovalue  =  []; 

index  =  index  -  1; 

indexvalue  =  1; 

%  Form  a  new  array 

for  m  =  1 : index 

for  n=  1 : Amplitudecount (m)  ; 

Histovalue (indexvalue)  =  Amplitudevalue (m) ; 
indexvalue  =  indexvalue  +1; 

end 

end 

Histolength  =  length (Histovalue) ; 

%  Recalculat  kurtosis  on  the  reduce  histogram  array 
ProdKurt  =  kurtosis (Histovalue) ; 


end 

end 

disp ( index) ; 
disp (ProdKurt) ; 
disp (X (index) ) ; 

%  Plot  RF  spectrum  and  the  threshold 
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figure  ( 3 ) ; 


plot (f /le6, a) ; 
axis  tight; 
grid; 

xlabel ( 'Frequency  (MHz)  ' )  ; 
ylabel ( 'Power  (dBm)  ' )  ; 

title  (strcat (filename (SELECT)  ,  --  ' , rbw,  ' 

Hz'),' Interpreter ' , ' none ' ) ; 

hold; 

xa  =  [f(l)/le6  f(end)/le6]; 
yb  =  [X (index)  X (index)]; 
plot (xa, yb, ' r-o ' ) ; 

%  Save  results  to  data  file 

hold  off; 

lab  =  num2str (SELECT) ; 

label  =  strcat (' Kurtosis-Figure ', lab, '. jpg ' ) 
saveas (gcf , label) ; 


end 
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Appendix  B.  Graphs  of  RF  Spectrum  Files  Calculated  Detection 

Threshold 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


ARL 

US  Army  Research  Laboratory 

CF 

crest  factor 

GHz 

gigahertz 

kHz 

kilohertz 

MHz 

megahertz 

PDF 

probability  density  function 

RBW 

resolution  bandwidth 

RF 

radio  frequency 
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